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We derive a phenomenological continuum saltation model for aeolian sand transport that can 
serve as an efficient tool for geomorphological applications. The coupled differential equations for 
the average density and velocity of sand in the saltation layer reproduce both known equilibrium 
relations for the sand flux and the time evolution of the sand flux as predicted by microscopic 
saltation models. The three phenomenological parameters of the model are a reference height for 
the grain-air interaction, an effective restitution coefficient for the grain-bed interaction, and a 
multiplication factor characterizing the chain reaction caused by the impacts leading to a typical 
time or length scale of the saturation transients. We determine the values of these parameters by 
comparing our model with wind tunnel measurements. Our main interest are out of equilibrium 
situations where saturation transients are important, for instance at phase boundaries (ground/sand) 
or under unsteady wind conditions. We point out that saturation transients are indispensable for a 
proper description of sand flux over structured terrain, by applying the model to the windward side 
of an isolated dune, thereby resolving recently reported discrepancies between field measurements 
and theoretical predictions. 



I. INTRODUCTION 

Aeolian sand transport, from the entrainment of sin- 
gle grains to the formation and movement of dunes, have 
been studied for a long time. One of the most impor- 
tant issues has been the relation q(u*) between the shear 
velocity u* and the saturated sand flux q. The simplest 
flux law, which gives a cubic relation between shear ve- 
locity and sand flux, was already introduced by Bagnold 
in 1936 |IJ]. Since that time, many new flux relations 
have been proposed and used by different authors. The 
most important improvement was to introduce a thresh- 
old to account for the fact that at low wind speeds no 
sand transport occurs. (A overview of the historical de- 
velopment can be found in Ref. One of the most 
widely used flux relations with threshold was proposed by 
Lettau and Lettau || . Analytical derivations of the flux 
relation starting from a microscopic picture deepened the 
understanding of the aeolian transport mechanisms a lot 
H-0. An application of sand flux relations are geomor- 
phological problems, where they are used to calculate the 
erosion rate from the wind shear stress in order to predict 
the evolution of a free sand surface or dune. However, all 
flux relations of the type g(tt*) assume that the sand flux 
is everywhere saturated. This condition is hardly fulfilled 
at the windward foot of an isolated dune, e. g. a barchan 
(crescent shaped dune, discussed in Section VII), where 
the bed changes rapidly from bedrock to sand. Corre- 
lated measurements of the sand flux and the wind speed 
performed by Wiggs || showed a large discrepancy be- 
tween the measured flux and theoretical predictions of 
the sand flux using the relation by Lettau and Lettau 
near the dune's windward foot. Numerical simulations of 
barchan dunes by Wippermann and Gross || that em- 
ploy this flux law also revealed this problem. Apart from 



the conditions at the dune's foot, it is conceivable that 
the sand flux may never reach saturation on the entire 
windward side of a dune, where the shear velocity in- 
creases gradually from the foot to the crest. Such effects 
are obviously not captured by an equilibrium flux law. 
To overcome the limitation of the equilibrium relations 
and to get information about the dynamics of the aeo- 
lian sand transport, numerical simulations based on the 
grain scale have been performed [jicj— . They showed 
that the typical time to reach the equilibrium state in 
saltation on a flat surface is approximately two seconds, 
which was later confirmed by wind tunnel measurements 
]l3[ . The problem of simulations on the basis of grains is 
that they can neither now nor in the near future be used 
to calculate the evolution of macroscopic geomorpholo- 
gies. 

In the following we derive a dynamic continuum model 
that allows for saturation transients and can thus be ap- 
plied to calculate efficiently the erosion in presence of 
phase boundaries and velocity gradients. In section || 
we introduce the phenomenology of aeolian sand trans- 
port. In section HI we develop a continuum model for a 



thin fluid like sand layer on an immobile bed including 
the time dependence of the sand transport and saturation 
transients. The following sections discuss special cases, 
where certain restrictions lead to simpler versions of the 
model. In section [V we discuss the saturated limit of 



the model and compare it with flux relations and experi- 
mental data from the litrature. In section^ we disregard 
the spatial dependence of the sand flux and concentrate 
on the time evolution of the saltation layer. In section VI 
we present a reduced "minimal model" that can easily be 
applied to geomorphological problems. Finally, we apply 
this model in section |VII| to predict the sand flux on the 
central slice of a barchan dune. 
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II. SAND TRANSPORT AND SALTATION 

Conventionally, according to the degree of detachment 
of the grains from the ground, different mechanisms of 
aeolian sand transport such as suspension and bed-load 
are distinguished. The bed-load can be further divided 
into saltation and reptation or creep. A detailed overview 
of this classification can be found in Ref. §. If we con- 
sider typical sand storms, when shear velocities are in the 
range of 0.18 to 0.6 ms -1 particles with a maximum 
diameter of 0.04-0.06 mm can be transported in suspen- 
sion. The grains of typical dune sand have a diameter 
of the order of 0.25 mm and are therefore transported 
via bed-load. For this reason we neglect suspension in 
the following discussion. Furthermore, we do not distin- 
guish between saltation and reptation but consider jump- 
ing grains with a mean trajectory length. 

The saltation transport can conceptually be divided 
into four sub-processes. To initiate saltation some grains 
have to be entrained directly by the air. This will be 
called direct aerodynamic entrainment. If there is al- 
ready a sufficiently large amount of grains in the air the 
direct aerodynamic entrainment is negligible and grains 
are mainly ejected by impacting grains. The entrained 
grains are accelerated by the wind along their trajectory 
mainly by the drag force before they impact onto the bed, 
again. This is called the splash process, which comprises 
the complicated interaction between bed and impacting 
grain and is currently the subject of theoretical and ex- 
perimental investigations |lj, fLS). Finally, the momen- 
tum transfered from the air to the grains leads in turn to 
a deceleration of the air. Through this feedback mecha- 
nism the saltation dynamics reaches an equilibrium state, 
characterized by a saturated density p s and an average 
velocity u s of the saltating grains. 

In the following we develop an effective continuum 
model in order to calculate the bed-load transport. The 
variables will be the density p and mean velocity u of 
the grains within a thin surface layer. The closed system 
of equations will have three phenomenological parame- 
ters. The first parameter a models the loss of energy 
in the splash process and characterizes the grain-bed in- 
teraction. It can be thought of as an effective restitution 
coefficient. The second parameter z\ is a reference height 
between the ground and the mean trajectory height and 
characterizes the air-grain interaction. These two pa- 
rameters determine the saturated sand flux q s , whereas 
the third parameter 7 determines the time scale T s or 
length scale l s of the saturation transients. 



III. A CONTINUUM MODEL FOR SALTATION 

We consider the bed-load as a thin fluid-like granu- 
lar layer on top of an immobile sand bed. This idea 
was introduced by Bouchaud et al. [|l6| and used in the 
following to model avalanches on inclined surfaces near 



the angle of repose ]l7, 18|. This general idea was later 
also applied to model the formation and propagation of 
sand ripples |l^,|2^]. To avoid cumbersome notations we 
restrict ourselves to a two dimensional description of a 
slice of a three dimensional system that is aligned with 
the wind direction. By this simplification we neglect the 
lateral transport, caused for instance by gravity or dif- 
fusion, which is typically an order of magnitude smaller 
than the flux in the wind direction. A further simplifica- 
tion jl6| [l8| is obtained by integrating over the vertical 
coordinate, which leads to scalar functions and equations 
for the moving layer. 

We start the derivation from the mass and momentum 
conservation in presence of erosion and external forces. 
Since the saltation layer can exchange grains with the 
bed, it represents an open system with the erosion rate 
r as a source term, 



d P , d 1 \ r 
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Here, p(x,t) and u(x,t) denote the density and velocity 
of the grains in the saltation layer, respectively. The ero- 
sion rate T(x,t) counts the number of grains per time 
and area that get mobilized. 

The most important forces acting on the grains are the 
drag force fdragix, t) when the grains are in the air, which 
accelerates the grains, and the friction force fbed(x, t) rep- 
resenting the complicated interaction with the bed, which 
decelerates the grains. Thus, we can write the momen- 
tum conservation for the saltation layer as 



du 



d_ 

dx 
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U = — (fdrag + fbed) ■ 



(2) 



Additional forces like gravity that might be important on 
inclined surfaces or diffusion caused by the splash process 
are neglected here to keep the model as simple as possi- 
ble, but could easily be added on the right hand side of 
Eq.(||). In the following sections we derive phenomeno- 
logical expressions for the erosion rate T and the forces 
fdrag and fbed, which will finally lead to a closed model. 



A. The Atmospheric Boundary Layer 

Sand transport takes place in the turbulent boundary 
layer of the atmosphere, near the surface. The Navier- 
Stokes equation /9air9tV + /5 Q i r (vV)v = — Vp+Vr, where 
v denotes the wind velocity, p a i r the density of the air, p 
the pressure, and r the shear stress of the air, reduces to 



dr 

dz 



or t — constant 



(3) 



by making the usual boundary layer approximation, ne- 
glecting d/dx,d/dy against d/dz and assuming steady 
d/dt = and horizontally uniform (vV)v = flow. At 
wind velocities for which sand transport is possible, the 
air flow is highly turbulent. Therefore, we can neglect the 
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bare viscosity of the air and identify r with the turbulent 
shear stress. The standard turbulent closure using the 
mixing length theory models the turbulent shear stress 
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where k w 0.4 denotes the von Karman constant. From 
equation (||) we obtain by introducing the characteristic 
shear velocity it», 



dv u* 
dz KZ 



with 



T 

Pair 



(5) 



Integration of Eq. (||) leads finally to the well known law 
of the wall for turbulent flow and therefore to the loga- 
rithmic profile of the atmospheric boundary layer, 

v(z) — — In — , (6) 

K Z 

where z denotes the roughness length of the surface. 



B. The Grain Born Shear Stress 

In presence of saltating grains near the ground, the 
air can transfer momentum to the grains, which thereby 
transport a part of the shear stress down to the surface. 
This idea was introduced by Owen |Ij and has widely 
been used in analytical saltation models and numerical 
simulations g-g|^|]J,|2l],|§. Accordingly, we distin- 
guish the grain born shear stress r g and the air born 
shear stress r a , which together have to maintain the over- 
all shear stress r, 



r = T a (z) + T g (z) = constant. 



(7) 



At the top z rn of the saltation layer, the air born shear 
stress T a has to be equal to the overall shear stress r, 

T a (z m ) = T. 

The typical trajectory of a saltating grain intersects 
with each height level z two times, once when ascending 
and once when descending. Between these two intersec- 
tions the wind accelerates the grain along its trajectory, 
thereby increasing its horizontal velocity between the as- 
cending and descending intersection. From this velocity 
difference we can calculate the shear stress transported 
by the grains at each level z, 

T g (z) = ®[u do wn(z) - Uup(z)} = ^Aug(z), (8) 

where $ denotes the flux of grains impacting onto the 
surface, and u up and Udown the horizontal velocity of the 
grains at the ascending and descending intersection of 
the trajectory with the height level z. 

As sketched in Fig. |l| we can relate both the horizon- 
tal sand flux q and the flux of grains impacting onto the 



surface $ to the number of grains N, their mass to, the 
saltation length I, and the flight T time for the trajectory, 



N m 
Taf 



q 

V 



(9) 



where a is the arbitrary width of the slice, considered 
here. Using Eq.(||) and (JsJ) the grain born shear stress 
can be expressed by the horizontal flux of grains, 



0) = -jAug(z). 



(10) 




FIG. 1. Sketch of the horizontal sand flux q caused by 
saltating grains passing the vertical rectangle. The dashed 
rectangle shows the surface area I times a that is used to cal- 
culate the flux <E> of grains impacting onto the surface, where 
I is the length of the saltation trajectory and T the flight time 
of the grain. 

In the following we restrict the discussion to the grain 
born shear stress r g o at the ground that is given by the 
momentum transfer from the grains to the bed during 
their impacts as sketched in Figure [j]. We denote quan- 
tities that are taken at the ground with an index 0, e.g. 
T gQ = Tg(zo), where zg is the height at the ground. We 
need the impact and ejection velocities of the grains or 
moreover the change in horizontal velocity Au 9 o at the 
ground to calculate the momentum transfer. In order to 
keep the discussion simple we will directly formulate the 
model in terms of mean values for the trajectory length I 
and its time T instead of writing everything in terms of 
not well known distribution functions. 

The average saltation length I and the average salta- 
tion time T are related by I = uT . We estimate both as 
the flight time and length of a simple ballistic trajectory, 



T = 



2u 



I 



2u z0 

t i 
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(11) 



where u z q is the vertical component of the initial veloc- 
ity of the grain and g the acceleration by gravity. The 
typical values for the time and length of a mean salta- 
tion trajectory depend on the shear velocity it*, T(tt*) f» 
1.7tt*/p and Z(it«) « I8v%/g fL4]| . For a shear velocity 
= 0.5 ms -1 we obtain T sa 0.08 s and I w 0.45 m. 
Calculations of grain trajectories confirm that this ap- 
proximation gives the correct order of magnitude for the 
flight time T, e.g. S0rensen obtained T = l.lhuJg [Ell. 
Inserting the saltation length of Eq.(|TT|) in Eq.(10) and 
using q — pu we obtain, 



T g0 = P 



9 AM g o 
2 u z0 



(12) 



3 



The relation between the impact and ejection ve- 
locity is in general defined by the splash function 

) giving the probability that a 
grain is ejected at a certain angle and velocity due to 
an imp acting grain with a certain angle and velocity 
0,|lJ^5|. The simplest possibility is to take the ver- 
tical component of the ejection velocity u z q proportional 
to the horizontal velocity difference Ak 9 o and to neglect 
the angle dependence, 



u z q = a Au 



gO- 



(13) 



Here, we have introduced the first model parameter a 
that represents an effective restitution coefficient for the 
grain-bed interaction. In principle, it can be calculated 
from the splash function S, but here we regard it as a 
phenomenological parameter to be determined later by 
comparing the model with experimental results. With 
Eq.(13|), Eq.(|l^) reduces to the simple result 
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(14) 



for the grain born shear stress. The fact that it is linear 
in the density p and independent of the velocity u is an 
interesting result that is not a priori obvious. Indeed, 
if we had chosen the effective restitution coefficient a (it) 
as a function of the velocity we would obtain a velocity 
dependence. Why is the velocity not or only of little im- 
portance for T g o? This can be explained by looking at 
the trajectory length of the saltating grains, which scales 
with the square of the horizontal velocity (I oc w 2 ) lead- 
ing to a decrease of the density of impacts with increasing 
u. This compensates the higher momentum transfer of a 
single grain due to its higher velocity. 



C. Erosion and Deposition Rates 

Grains impacting onto a bed of grains either rebounce 
directly or remain on the bed. In the latter case the en- 
ergy of the impact can dislodge several new grains lead- 
ing to a strong amplification of the grain density in the 
saltation layer and in turn to erosion. The probability 
that a certain number of grains leaves the surface with a 
certain velocity is in general given by the splash function 
S (nj |l4|, [HI , mentioned above. A part of the informa- 
tion contained in the splash function has already been 
comprised into the model parameter a that relates the 
average horizontal velocity difference to the vertical com- 
ponent of the average ejection velocity. Another part can 
be represented by the average number n of grains disloged 
by an impacting grain. Together, a and n determine the 
simple effective splash function used in this model, 



S(u z0 , Au g0 ) = nS(aAu g0 - u z0 ). 



(15) 



Using the average number n of grains dislodged by an 
impacting grain we can define the erosion rate as the net 



average rate of grains leaving the surface, which is the 
difference between the flux of grains leaving the surface 
and the flux of grains impacting onto it, 



r = 1) = 



T ffO 

AUgO 



(n-1). 



(16) 



It is important to realize that the air born shear stress 
within the saltation layer, and therewith n, is lowered if 
the number of grains in the saltation layer increases and 
vice versa. This is the feedback effect discussed in sec- 
tion E3, According to Owen , the air born shear stress 
at the bed r a o in the equilibrium is just large enough 
to keep saltation alive and therefore close to the thresh- 
old r t = Pairul t . For T a Q > T t , the number of grains in 
the saltation layer increases in a chain reaction [n > 1) 
whereas for t q o < r t (n < 1) saltation cannot be main- 
tained. To model the average number n of dislodged 
grains out of equilibrium we write n as a function n(r Q /r t ) 
with n(l) — 1. We furthermore assume that n can be ex- 
panded into a Taylor series at the threshold and neglect 
all terms after the linear order, 



TgQ 



= 1 



, Ta0 1 



(17) 



The model parameter 7 characterizes the strength of the 
erosion and determines how fast the system reaches the 
equilibrium or reacts to perturbations. It depends on 
microscopic quantities such as the time of a saltation 
trajectory or the grain-bed interaction, which are not 
available in the scope of this model. Therefore, we have 
to determine 7 later by comparison with measurements 
or microscopic computer simulations. 

If we insert Eq.(^) and ( p"7j ) in Eq.(16), we obtain for 
the erosion rate 



AUgQ 



TgQ 



- 1 



(18) 



Assuming that the difference between impact and eject 
velocity of the grains is proportional to the mean grain 
velocity (Au g o oc u) finally leads to 



TgQ 



TgQ 



Tt 



(19) 



where the proportionality constant is incorporated in 7. 
In order to close the system of equations we have to in- 
sert the expression for the grain born shear stress r g o 
from Eq.(|lJ). 

Up to now, we discussed the entrainment of grains due 
to impacts of other grains, but if the air shear stress 
exceeds a certain value r ta > r t , called aerodynamic en- 
trainment threshold, grains can directly be lifted from 
the bed. These directly entrained grains have been ne- 
glected up to now, because they are only important to ini- 
tiate saltation p0| , pTlj • Anderson Q proposed an aero- 
dynamic entrainment rate proportional to the difference 
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between the air born shear stress r a and the threshold is now varying in z, 



Tta, 



Tta 



Tg 



T/a 



(20) 



where $ a w 5.7 10 _4 kgm _2 s _1 [0 is a model param- 
eter defining the strength of the erosion rate for aero- 
dynamic entrainment. This formula for the direct aero- 
dynamic entrainment rate T a has a similar structer as 
Eq.(|l9|) for saltation induced entrainment, but the pref- 
actors arc different. 



D. Forces 



In section [II we introduced the drag and friction 
forces, fdrag and fbed, acting on the saltation layer. In 
the following we have to specify these forces. Modeling 
the friction force is simple, because we already derived an 
expression for the grain born shear stress at the ground 
T g0 , Eq.(Q[). The bed friction ff, e d has exactly to com- 
pensate this grain born shear stress, 



fbed — — t. 



'go- 



(21) 



We represent the wind force acting on the grains inside 
the saltation layer by the Newton drag force, 



1 



ird 2 



Fdrag — Paired . \pair ^g)\^a 



(22) 



of a spherical particle, where d denotes the grain diam- 
eter, Cd the drag coefficient, v g the velocity of a grain, 
and v a i r the velocity of the air. Multiplying Fd ra g with 
the density p of the saltation layer and dividing it by the 
mass m of a grain with diameter d and density p qU artz, 
we obtain the drag force acting on a volume element of 
the saltation layer, 



r 3 p a i r 1 . 

J drag = Pj^d~ "j {V e ff 

^ Pquartz W 



' J eff 



(23) 



Here, v e ff is an effective wind velocity, which is the wind 
speed taken at a reference height Z\ within the saltation 
layer. This reference height is another model parameter 
and has to be determined by comparing the sand flux 
with measured data as we will do it in section IV. Ne- 



glecting the effect of the saltating grains on the wind 
field, we could use the logarithmic profile at z\ to calcu- 
late the effective wind speed. But saltating grains in the 
air change the air shear stress and the wind speed. This 
is the feed back effect mentioned above and the mecha- 
nism how an equilibrium sand flux is reached. In order 
to calculate the perturbed wind speed we use again the 
standard turbulent closure ([s]), which relates the strain 
rate to the turbulent shear stress. In contrast to the case 
without grains, where the shear stress r of the air is con- 
stant in z, the air born shear stress r a as given by Eq.(|^) 



dv 
dz 



KZ 



T g( z ) 



(24) 



The profile of t (z) was found to be nearly exponential 
in simulations [p.l| . Furthermore, we already know the 
grain born shear stress at the ground r g o- Hence, we can 
write 



dv 

dz 



KZ 



1 



TgO 



I* 



(25) 



where z m denotes the mean saltation height. The inte- 
gration of Eq.(|25|) has to be performed from the rough- 
ness height zq to a reference height z\ < z m . Therefore 
we linearize the exponential function to integrate Eq. (|2E 
and obtain the effective wind velocity 



V eff 



with 



K 
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(26) 



(27) 



For a reference height z\ much smaller than the mean 
saltation height z m (zq < z\ <ti z m ) we can simplify 
Eq.© to 



' J eff 



K 



TgO 

T 
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(28) 



For vanishing grain born shear stress T g o — > 0, the ef- 
fective wind velocity v e f f reduces to the velocity of the 
undisturbed logarithmic profile at height z\. The values 
of the parameters zo and z m can be obtained from mea- 
surements, whereas the value of the reference height z\ 
is a free phenomenological parameter of the model that 
we have to estimate by comparing the saturated flux pre- 
dicted by our model with measurements. 



E. Closed Model 

So far, we have introduced the erosion rate T, the 
drag force fdrag, and the bed interaction f^d and ex- 
pressed them in the preceding sections in terms of the 
mean density p and the mean velocity u of the saltating 
grains. Furthermore, we have introduced two model pa- 
rameters a and z\ determining the equilibrium state of 
the saltation layer, and the parameter 7 that controls the 
relaxation to equilibrium. Inserting Eqs.(|l9|) and ( [Til ) in 
Eq.(|l]) leads to an equation for the sand density p in the 
saltation layer, 



dp d 



7 



2a r t 



Tt - P -\l- 



1 



2a t — T t 
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(29) 
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Here, we can identify two important physical quantities, 
the saturated density p s and the characteristic time T s 
that define the steady state and the transients of the sand 
density p, respectively, 



Pa 



2a 



T„ = 



2au 



9 l{r-r t )' 



(30) 



(31) 



With these expressions we can rewrite Eq. ( |2^ ) in a more 
compact form, 



dp d If p 



(32) 



Direct aerodynamic entrainment, i.e. the initiation of 



saltation, has been neglected in Eq.(32), but can easily 
be included by adding the erosion rate T a of E q.(p0| ) to 
the right hand side. Furthermore, inserting Eqs.(|21|) and 
( p3| ) in Eq. (||) leads to a model for the sand velocity u in 
the saltation layer, 



du 
dt 



d_ 

dx 



« ^ pair i 
A Cd ~n H^ eff 



u)\v eff 
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(33) 



with vefj defined in (|28|). Finally, the Eqs.(|32j), (|30|), 
(|3l|), (j33[), and ( p8| ) define the closed model for the sand 
flux in the saltation layer. 

We want to emphasize that T s (r, u) and l s (r, u) — T s u 
are not constant, but depend on the external shear stress 
r of the wind and the mean grain velocity u. Using 
Eq.([i"T|) we can relate the characteristic time T s and 
length l s of the saturation transients to the saltation time 
T and the saltation length I of the average grain trajec- 
tory, 



T, =T- 



n 



l(r-n) 



L = l- 



Ti 



7(T - T t ) 



For typical wind speeds, the time to reach saturation is in 
the order of 2 s [p~0| — p~2|] . Assuming a grain velocity of 3- 
5 m s _1 we obtain a length scale of the order of 10 m 
for saturation. This length scale is large enough to play 
an important role in dune formation. We want to em- 
phasize that the characteristic length l s of the saturation 
transients, Eq.(^), naturally result from the saltation 
kinetics (in contrast to the heuristic "adaptation length" 
of Ref. [p4fl ). Their functional dependence can be in- 
terpreted in the following way. The dominant mecha- 
nism to adapt the grain density in the saltation layer 
due to a change in external conditions is by the chain 
reaction process modeled in Eq.(|l7j). Hence, l s and T s 
depend inversely onto the "stiffness" j(r/Tt — 1) of this 
multiplication process. Since the average grain can only 



influence the number of grains in the saltation layer at 
the discrete times and positions where it interacts with 
the bed, the saturation time/length are proportional to 
the time/length of the average trajectory. The resulting 
non-trivial dependence of the saturation kinetics on the 
external conditions, which may be appreciated from Fig- 
ure ^| in Section VI , is essential for a proper description 
of aeolian sand transport in structured terrain. 

The mass and momentum Eqs.(|32|) and ( |33] ) are cou- 
pled partial differential equations and difficult to solve. 
In the following sections we will first discuss the fully 
saturated situation and later some dynamical properties 
of the saltation layer. Finally, we simplify the model in 
order to arrive at a minimal model that can easily be 
applied to geomorphological problems. 



IV. SATURATED FLUX 

The full dynamics of the saltation layer must be eval- 
uated numerically, whereas the saturated flux — the sta- 
tionary solution (d/dt = 0) for a constant external shear 
stress (r(x,t) = r) and a homogeneous bed(d/dx = 0) 
— can be calculated analytically from Eq.(p9|) and (|3^). 
For shear velocities below the threshold (it* < it* t ) the 
solution is trivial, 

p s (u*) — 0, u s (m*) = 0, and q s {u*) = 0. (35) 

Above the threshold (tt* > u*i) we obtain from Eq.(|30|) 
for the steady state density p s , 



2ap a 
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(36) 



Likewise we obtain from Eqs.(|33|) the steady state veloc- 
ity u s , 



u s (u*) 




" ' ' C 2 U *t (0T\ 

— o" \-Ust, (37) 



(34) where 



u st = u s (u*t) = m 

K Z 



>2gdp qv 



3aCd Pai 



(38) 



is the minimum velocity of the grains in the saltation 
layer occurring at the threshold. In contrast to the den- 
sity p s the velocity u s does not go continuously to zero 
near the threshold it*t. This is intuitively obvious, be- 
cause grains at the threshold have already a finite veloc- 
ity u s t- Finally, we can write for the steady state flux 
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Ust 



(39) 
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For large wind speeds (it* 3> lt*j) the flux is asymptot- 
ically proportional to u 3 ,, which is in accord with the 
predictions of other saltation models [p|,|2f|, 25 j. 

The saturated flux q s , given by Eq.(j39|), is now used 
to determine the model parameters a and Z\ by fitting 
it to flux data measured in a wind tunnel by White 
p6fl . Using the literature values: g — 9.81 ms~ 2 , p air = 
1.225kgm" 3 , p quar t z = 2650kgm~ 3 , z m = 0.04m, 
z = 2.5 10~ 5 m D = d = 250 pm, C d = 3 and 
u*t = 0.28 ms -1 |g,||,[Ll| we obtain for the two model 
parameters a = 0.35 and Z\ = 0.005 m. For comparison 
we fitted to the same set of data the sand transport laws 
given by Bagnold [B5| , 



„ Pair Id 3 
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Lettau and Lettau |p|, 
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(40) 



(41) 



and S0rensen pi], 



„ Pair , 
9S = (~>S M*(M„ 
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u*t)(u* +7.6*w* t + 205) (42) 



and obtained C B = 1-98, C L = 4.10, and C s = 0.011. 

The results are shown in Figure |^, where the flux nor- 
malized by qo = Pair/igul) is plotted versus the nor- 
malized shear velocity u*/u*f Our result resembles 
S0rensen's equation but differs from the flux relation 
given by Lettau and Lettau. For high shear velocities 
all transport laws show a cubic dependence on the shear 
velocity. 
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FIG. 2. Comparison of the different theoretical flux rela- 
tions (|3|)-@ fitted to wind tunnel data of White. The fluxes 
are normalized by qo = Pair/ (<?«»)■ The saturated flux of our 
model and the relation of S0rensen reproduce quite well the 
data, whereas the others do not show the same structure. 



The comparison of the saturated sand flux with ex- 
perimental data determined the two phenomenological 
parameters a and z\. As anticipated above, the value ob- 
tained for z\, the reference height for momentum trans- 
fer, is well below z m , the mean height of the saltation 
trajectories. This justifies a posteriori the linearization 
in Eq.(|2F 



V. DYNAMICS 

After the saturated case has been studied in the pre- 
ceding section, we now investigate the dynamics of the 
saltation layer in order to get an estimate for the sat- 
uration time T s and thus for the model parameter 7. 
Figure || shows numerical solutions of Eq.(p2[) and d33| ) 
for the time evolution of the sand flux q = u p using spa- 
tially homogeneous conditions (d/dx = 0). To get rid of 
the free parameters in Eq.(p0|) we neglected r a , thus dis- 
regarding direct aerodynamic entrainment, and assumed 
instead a small initial density. Due to the multiplication 
effect of the saltation process, the flux increases first ex- 
ponentially and reaches the equilibrium state after pass- 
ing through a maximum at t w 2 s. The time transients 
are controlled by the parameter 7 and compare well with 
measurements by Buttcrfield ]13j and microscopic simu- 
lations by Anderson @|ll| and McEwan flf] for 7 w 0.4. 
An important result of the simulations of Anderson was 
the dependence of the saturation time on the shear ve- 
locity and the overshoot near t w 2 s. Both features are 
well reproduced by our model. 
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FIG. 3. Numerical simulations of the time evolution of the 
full model given by Eq.(^2|) and ( pS3| ) with a constant shear 
velocity u*. 



VI. A MINIMAL MODEL FOR 
GEOMORPHOLOGICAL APPLICATIONS 

The change of desert topographies, e.g. the movement, 
growth, and shrinkage of dunes, depends mainly on the 



7 



sand transport or erosion and the perturbations of the 
wind field caused by the topographies themselves. Here, 
we restrict ourselves to the problem of the sand transport 
and assume that the shear velocity above a certain topog- 
raphy is known. The time evolution of the topography 
h(x, t) is then given by the mass conservation, 



dh 

at 



Ps 



1 dq 

dx" 



(43) 



where p san d is the mean density of the immobile dune 
sand. To obtain the sand flux q(x, t) one can in princi- 
ple solve the coupled differential equations for the den- 
sity p, Eq. (jj2|) , and velocity u, Eq.(|3^), of the sand in 
the saltation layer. However, for most geomorphologi- 
cal applications a simplified version of our model will be 
sufficient. In the following, we first deri ve th is "mini- 
mal model" from the Eqs. given in Section HIE and then 
show its usefulness b y dis cussing a particular practical 
application in Section VII. 

The first simplification is to use the stationary solution 
(d/dt = 0) of Eq.© and @. This can be justified by 
the fact that there are several orders of magnitude be- 
tween the time scale of saltation (approximately 2s) and 
the time scale of the surface evolution of a dune (several 
days or weeks). 

Next, we consider the convective term (ud x u) that is 
only important at places where large velocity gradients 
occur. This is for instance the case in the wake region 
behind the brink of a dune, where the wind speed at the 
ground decreases drastically due to the flow-separation 
of the air. Here, the inertia of the grains becomes im- 
portant. To solve the model for the deposition in a wake 
region we want to consider an idealized brink situation, 
where both the wind speed and the friction with the bed 
drop discontinuously from a finite value to zero. In this 
case, Eq.(p3[) reduces to 



Idu 2 
2~dx 



3 „ Pair 1 2 

A Cd 1 U ■ 

P quartz & 



(44) 



This predicts an exponential decrease of the grain 
velocity over a characteristic length scale Idep = 
4<i p quartz I '(3Cd Pair) ~ 0.25 m. Hence, the deposition 
takes place within a length Idep / s much shorter than 
the saturation length l s on the windward side. Field mea- 
surements of lee side deposition agree with this conclu- 
sion ^7|. Outside the wake regions, on the other hand, 
we can neglect the convective term (ud x u) and we obtain 
the mean stationary grain velocity from Eq. (j33|) 



U (p) = V eff(p) 



*2gdp qv 



3aCdp a 



(45) 



where u e //(p) is defined by Eq.(p8|) and (pT[). Further- 
more, near the threshold r t we can approximate v e ff(p) 
by v e f f {p s ) making only a negligible error. For high shear 
stresses this is not in general possible. But, the sand flux 



in macroscopic geomorphological applications is nearly 
everywhere saturated, apart from places where external 
variables change discontinously, e.g. at a phase bound- 
ary bedrock/sand or at a flow-separation (see Fig. [?]). 
Hence, we can replace the density p by the saturated den- 
sity p s for most applications, where the exact value of the 
flux at these places is not of importance. The advantage 
of this approximation is that the velocity u is decoupled 
from the density p and we can insert the saturated ve- 
locity u s of Eq.© in Eq.©. Rewriting Eq.@ in in 
terms of the sand flux q = pu s and the saturated sand 
flux q s — p s u s of Eq. (|3^) leads to an equaation for the 
sand flux q, 



d 1 A q 
Tp9 = T q l 

Ox l s \ q s 



(46) 



where 



2a ui 



n 



7 9 T-T t 



(47) 



is the saturation length depicted in Figure |5|. 

We want to emphasize that our most im por tant re- 
sult for geomorphological applications is Eq.poj), which 
extends a common saturated flux law by incorporating 
saturation transients. It is to some extent independent 
of the particular form of the functions / s and q s given in 
Eq.([47|) and (|39|). The latter can be regarded as addi- 
tional predictions of our model which one could also re- 
place by other phenomenological relations or data tables 
obtained from wind tunnel measurements if this turns 
out to be more suitable. 



0.09 
0.08 
0.07 
0.06 
0.05 



M 0.04 



0.03 
0.02 
0.01 



= 0.3 
= 0.4 
= 0.5 
= 0.6 
= 0.7 



01 23456789 10 

x in m 

FIG. 4. Numerical solution of the sand flux equation 
for different shear velocities u*. The model parameter 7 = 0.2 
that defines the length and time of the saturation transients 
was chosen here so that saturation is reached between 1 s and 
2 s. Due to the simplifications made with respect to the full 
model, the initial over-shoot is lost. 
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FIG. 5. The saturation length l s from Eq.(^) is asymptot- 
ically constant for high shear velocities (horizontal line), but 
diverges for shear velocities near the threshold. 

The saturated flux q s ( |3^ ) and the model parameters 
a = 0.35 and z\ = 0.005 m have been discussed in sec- 
tion IV and can be used unchanged in Eq.(|6|). But with 



respect to the full model of section [v| the value of 7 gets 
renormalized due to the simplifications made in order to 
obtain Eq. ( E6|) . The most obvious difference of the so- 
lution of Eq. (|46|) , depicted in Figure [|, compared to the 
result of the full model in section ^ is the missing over- 
shoot. The value of 7 = 0.2 had to be adapted in order 
to obtain saturation transients between 1 s and 2 s for 
typical values of u*. The fact that the saturation length 
l s increases for shear velocities near the threshold is un- 
changed and shown in Figure ^ and 0. 



VII. 



THE SAND FLUX ON THE WINDWARD 
SIDE OF A DUNE 



windward side 



brink 




• horns 

FIG. 6. Sketch of a barchan dune. 
The mass conservation, Eq.(p3[), and the saturated 



sand flux relation by Lettau and Lettau, Eq.(41), have 
been used many times to predict the evolution of a dune 
|^,|2^,|9). The most successful work was done by Wip- 
pcrmann and Gross P], but even they obtained an un- 
physical deposition at the windward foot of the dune, 



which tended to stretch the dune in length and to flat- 
ten it, finally. This was avoided by an ad hoc smooth- 
ing operation, which led to a numerically stable sim- 
ulation by hiding the problem at the phase boundary 
(bedrock/sand). The boundary problem is evident for 
an isolated dune like a barchan sketched in Figure ||. Re- 
cent field measurements || of the wind speed and the 
sand flux on the central slice of a barchan dune show a 
large discrepancy between the measured sand flux at the 
windward foot and the sand flux predicted by Eq.(^l[) for 
the measured wind speed. The measured sand flux was 
a monotonously increasing function, whereas the calcu- 
lated sand fluxes decreased near the dune's foot due to 
the depression of the wind velocity. A decrease of the flux 
is correlated with deposition at the dune's foot leading to 
a flattening of the dune. Already from the measurements 
it is evident, that a saturated flux law cannot be applied 
near a phase boundary, where the bed changes suddenly 
from bedrock to sand. A further problem apears if we in- 
tegrate Eq.(p3|) in time using a flux relation of the form 
q(u*) that determines the sand influx qi n = q[u^(xi n )] 
and outflux q out = q[u*(x ou t)] at the boundaries. This 
might be true for the outflux, but not for the influx, which 
depends normally on the upwind conditions and not on 
the wind speed at the boundary. In particular, the influx 
of an isolated dune in a dune field should depend on the 
outflux of several dunes upwind. Even situations without 
influx are possible, e.g. if there is vegetation around the 
dune. 

To elucidate this problem further, we calculated the 
shear stress using FLUENT 5 |}0| with a /ce-turbulence 
model for the central profile, parallel to the wind direc- 
tion, of a barchan measured in Morocco. (For more de- 
tails see dune no. 7 in [plf.) The exact calculation of the 
flow field is not of importance for the following discus- 
sion, which only relies on characteristic qualitative fea- 
tures such as the depression near the dune's foot. For 
the calculation of the sand flux we assume bedrock up 
to the dune's foot (xi„=25m), where erosion is not pos- 
sible (g=constant). The sand flux over the bedrock is 
therefore equal to the influx at the boundary Xi n . The 
sand flux q^ according to the relation by Lettau and 
Lettau, Eq.(f4l|), and the solution of the saturated flux 
q s , Eq. (p9[) , predicted by our model are depicted in Fig- 
ure [f| Both models exhibit the unphysical deposition 
at the dune's foot described above. This problem is re- 
solved, when the "minimal model", Eqs.(p6|), ([47j), and 
(|39|), is applied. In contrast to a saturated flux law, the 
boundary conditions can freely be chosen in this model. 
Here, we used a constant influx q.i n , much smaller than 
the saturated one, which represents the inter-dune sand 
flux. At the outflow boundary we applied dq/dx = 0. 
The solution is plotted in Figure (7] in comparison with 
the predictions of the simple flux relations. Due to the 
saturation transients we obtain a monotonously increas- 
ing flux and therefore no deposition of sand at the dune's 
foot. Finally, we want to point out that the flux on the 







entire windward side is never fully saturated due to the 
continuously increasing shear stress. However, away from 
the boundary it is always close to the saturated value 
(compare q and q s in Fig. (?]), which justifies the simpli- 
fications made in the "minimal model" . 
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FIG. 7. Top: shear stress r of the air at the surface, cal- 
culated by a fce-turbulence model using FLUENT 5. Center: 
The sand flux according to the "minimal model", Eq.(|46|), the 
saturated sand flux q s , Eq.(|3S|), and the sand flux qL predicted 
by the Lettau and Lettau relation ( ftl| ) have been calculated 
using the shear stress depicted at the top. Bottom: Height 
profile of the symmetry-plane of a barchan. 



VIII. CONCLUSION AND OUTLOOK 

We derived a phenomenological sand transport model 
that reproduces the equilibrium sand fluxes measured in 
wind tunnels. The predicted evolution of the sand flux 
with time has the same qualitative average behavior as 
the sand flux calculated with saltation models working 
on the microscopic grain scale. Finally, we proposed a 
"minimal model" that can be used as an efficient tool for 
geomorphological applications, such as the formationn 
and migration of dunes. Furthermore, we showed that 
this model extends in a general way common saturated 
sand flux relation to a model that incorporates saturation 
transients. 

The phenomenological parameters for the saltation 
model, which are an effective restitution coefficient, a 



reference height within the saltation layer, and a satu- 
ration length have been estimated by comparison with 
measurements. 

Using these parameters we applied our model to a 
geomorphological problem and calculated the sand flux 
over a dune. We emphasized the importance of a non- 
equilibrium flux relation for the correct modeling of phase 
boundaries, e.g. bedrock/sand, as they naturally occur 
for isolated dunes. Furthermore, we have shown that the 
model predicts saturation transients to persist over the 
entire windward side of a dune, where the shear stress 
increases from the foot to the brink. Therefore, we claim 
that the saturation length defines a length scale that is 
important for dune morphology in general. The inves- 
tigation of the implications of saturation transients for 
dune formation will given in Rcf. |32j. In particular, 
the question of shape differences between small and large 
dunes or the minimum size for slip-face formation are of 
interest and are discussed there. The extension of the 
model to three dimensions and inclined surfaces, which 
is necessary for applications to arbitrary terrains will be 
the subject of futur work. 
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